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NUMERICAL  TREATMENT  OF  DIFFERENTIAL  AND  INTEGRAL 
EQUATIONS  BY  THE  p  AND  h-p  VERSIONS  OF  THE  FINITE 

ELEMENT  METHOD 

(Continuation  of  ‘Analysis  of  the  Performance  of  Mixed  Finite  Element  Methods’) 

Principal  Investigator:  Manil  Suri 
Co-investigator:  Christoph  Schwab 
Grant  AFOSR  89-0252 

Final  Scientific  Report-January,  1992 
Abstract 

Our  research  can  be  divided  into  three  main  groups. 

1.  An  investigation  into  the  problem  of  “locking”  which  arises  in  the  approxima¬ 
tion  of  parameter-dependent  problems.  We  have  developed  a  general  theoretical 
framework  to  analyze  this  phenomenon  and  have  characterized  the  locking  and 
robustness  of  different  finite  element  schemes  for  various  problems. 

2.  Results  on  the  p  and  h  —  p  versions  of  the  finite  clement  method. 

.  (a)  Optimal  approximation  results  for  the  p  version  BEM  in  three  dimensions. 

(b)  Analysis  of  a  p  version  mixed  method  for  quasilinear  problems. 

(c)  Analysis  of  the  error  in  Z,2  and  lower  norms. 

(d)  Investigation  of  quadrature  schemes  and  related  errors. 

3.  Results  obtained  by  C.  Schwab. 

(a)  Singularities  of  solutions  for  the  equations  of  three  dimensional  elasticity 
and  hydrodynamics  in  domains  with  edges  and  vertices. 

(b)  Numerical  evaluation  of  singular  surface  integrals  in  the  boundary  element 
method. 

(c)  Calculation  of  optimal  shear  correction  factors  for  plate  models. 
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1.  Introduction 


The  goals  of  this  project,  as  described  in  the  original  proposal  which  was  funded,  were 
to  investigate  the  p  and  h—p  versions  of  the  finite  element  method  in  various  areas,  such 
as  the  numerical  approximation  of  three-dimensional  PDEs  ana  integral  equations,  the 
investigation  of  mixed  methods  for  these  versions  and,  most  importantly,  the  uniform 
approximation  of  parameter-dependent  problems  by  these  versions.  By  the  p  version, 
we  mean  a  finite  element  method  where  the  mesh  is  kept  fixed  and  the  polynomial 
degrees  are  raised  to  increase  accuracy.  In  the  h  —  p  version,  both  the  mesh  and  the 
polynomial  degrees  are  changed.  Research  on  these  methods  (which  has  received  strong 
support  from  AFOSR,  among  others)  has  led  to  various  important  results  which  show 
that  the  p  and  h  —  p  versions  are  particularly  suitable  for  linear  elliptic  PDEs  (like 
those  arising  in  structural  analysis),  with  extremely  high  rates  of  convergence  being 
obtained.  Various  commercial  codes  (used  in  the  design  of  aircraft,  for  example)  have 
been  developed  on  the  basis  of  such  research  -  we  mention  PROBE  (developed  originally 
by  Noetic  Technologies,  MO  and  now  part  of  MSC  software)  and  STRIPE  (developed 
by  the  Aeronautical  Research  Institute  of  Sweden). 

The'  first,  and  most  important  part  of  our  research  has  been  on  the  approximation 
of  parameter-dependent  problems.  This  can  often  suffer  from  a  phenomenon  called 
“locking”,  which  causes  the  deterioration  of  certain  finite  element  approximations  for 
some  ranges  of  the  parameter.  Locking  can  cause  computational  results  in  real-life 
situations  to  turn  out  to  be  completely  inaccurate  and  has  been  described  by  Richard  H. 
MacNcal  (chairman,  MacNcal-Schwcndlcr  Corporation),  as  “the  most  virulent  problem 
that  afflicts  lower-order  schemes”  (ICIAM  91  lecture).  Our  goal  has  been  to  set  up 
a  mathematical  framework  to  define  and  analyze  locking,  our  suggested  strategy  to 
overcome  its  effects  has  been,  in  general,  to  use  higher-order  schemes.  In  this  context, 
we  have  shown  that  the  p  and  h  —  p  versions  can  be  particularly  effective,  as  can 
be  higher-order  h  version  schemes.  This  research  is  a  starting  point  for  many  of  the 
problems  we  will  be  doing  as  part  of  our  new  grant. 

Our  second  area  of  research  deals  with  results  obtained  for  the  p  and  h  —  p  versions 
of  the  finite  element  method  on  the  other  topics  mentioned  above  and  some  additional 


related  topics.  We  have  analyzed  a  mixed  p  version  finite  clement  method  for  a  quasi- 
linear  elliptic  problem,  obtaining  optimal  results.  We  have  investigated  the  optimality 
of  the  p  version  in  lower  order  norms,  when  applied  to  domains  with  corners.  We 
have  also  obtained  several  results  on  the  accuracy  of  the  p  version  in  the  presence  of 
numerical  quadrature,  both  for  source  and  eigenvalue  problems.  Since  January,  1991, 
Dr.  Christoph  Schwab  has  also  been  funded' by  this  grant.  In  joint  research,  we  have 
investigated  the  optimal  rate  of  convergence  of  the  p- version  boundary  element  method 
(BEM)  over  polyhedral  surfaces  in  R 3  and  obtained  sharp  rates  of  convergence,  the  first 
results  of  this  kind  in  three  dimensions. 
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In  addition,  Dr.  Schwab  has  investigated  jointly  with  W.L.  Wendland  (Germany) 
the  numerical  quadrature  of  singular  and  hypersingular  integrals  over  curved  surfaces, 
which  arise  for  example  in  boundary  element  methods  for  three  dimensional  crack 
problems.  In  other  work,  using  the  Fourier  analysis  of  hierarchical  plate  models,  a 
mathematically  rigorous  framework  for  the  optimal  selection  of  the  so-called  shear 
correction  factors  in  the  Reissner  Mindlin  and  related  plate  models  has  been  given. 
Finally,  in  joint  work  with  V.A.  Kozlov  and  V.G.  Maz’ya  (Sweden),  the  vertex  singular 
exponents  for  the  Lame-  and  Stokes-  equations  in  three-  dimensional  domains  with 
vertices  and  edges  have  been  analyzed  and  sharp  estimates  in  dependence  on  the  local 
vertex  geometry  obtained. 


2.  Results  achieved 

In  this  section,  we  describe  the  results  achieved  under  the  grant  over  the  past  three 
years. 

2.1.  The  problem  of  locking  in  the  FEM 

The  numerical  approximation  of  various  parameter-dependent  problems  may  suffer 
when  the  parameter  t  lies  close  to  a  limiting  value,  a  phenomenon  called  locking.  For 
example,  consider  the  case  of  an  elastic  plate.  As  the  thickness  tends  to  zero,  the  so¬ 
lution  Of  the  3-d  elasticity  problem  tends  to  that  of  the  Kirchoff  plate  and  Kirchoff’s 
hypothesis  gets  imposed  as  a  constraint  on  both  the  exact  and  on  the  approximate  solu¬ 
tion.  Locking  occurs  when  there  are  not  enough  functions  in  the  approximate  subspace 
satisfying  this  constraint  to  maintain  the  level  of  approximability  expected.  A  similar 
efTcct  occurs  when  the  Poisson  ratio  v  is  close  to  0.5,  in  which  case,  the  constraint  of 
incompressibility  can  cause  problems. 

Locking  effects  can  have  extreme  consequences  in  terms  of  “real-life”  computations. 
There  have  been  many  engineering  and  mathematical  papers  which  describe  various 
locking  effects  and  analyze  methods  (mainly  mixed  methods)  to  overcome  it.  However, 
the  treatment  :n  these  papers  has  been  ad  hoc,  with  the  term  “locking”  being  used  to 
describe  very  different  phenomena  (compare  (4]  with  [5],  for  example). 

In  [P5],  we  have  developed  a  systematic  theory  of  locking  which  includes  mathemat¬ 
ically  precise  definitions  and  also  both  asymptotic  and  practical  methods  of  quantifying 
the  locking  and  robustness  of  various  schemes.  Our  theory  provides  a  general  frame¬ 
work  for  the  analysis  of  the  examples  listed 'in  the  above  references.  An  important 
subclass  of  problems  we  consider  is  one  where  the  parameter  t  enters  the  variational 
problem  in  the  following  way:  find  ut  £  V  satisfying 

a(ut,v) +^(Cut,Cv)  =  Ft(v)  Vv  €  V. 


(2.1) 


4 


Here,  a(>,  •)  is  coercive  on  V  and  C  is  a  bounded  operator  acting  on  V.  As  t  — *■  0, 
the  solution  ut  satisfies  the  constraint 


(2.2) 


Cuo  =  0. 


For  example,  if  (2.1)  represents  the  2-d  elasticity  equations  and  t  »  1  —  2t/,  where  v 
is  the  Poisson  ratio,  then  Cu  =  divu  and  the  limiting  case  of  (2.1)  when  v  =  1/2 
is  Stokes’  problem,  whose  solution  is  divergence  free.  Similarly,  in  the  case  of  the 
Reissner-Mindlin  plate, 

Cu  =  C(<f>,uj)  =<f>  —grad  a; 

(where  (f>  measures  the  shift  in  the  normal  and  u>  the  transverse  displacement),  so  that 

(2.2) ,  as  mentioned  above,  is  simply  KirchhofF’s  constraint. 

For  problems  of  the  form  (2.1),  we  may  often  decompose  the  solution  as 

(2.3)  ut  =  xi0  +  (u<  -  u0)  =  «o  +  z< 


where  «o  is  the  solution  of  the  limit  problem  (t=0).  A  necessary  condition  for  the 
absence  of  locking  is  that  u0  be  optimally  approximated  under  the  constraint  (2.2),  i.e. 


(2.4)  inf J|u0  —  vllv  «  inf  ||u0  ~v\\v 

v£SN  ,esN 

. .  Cv=0 

where  SN  C  V  is  the  approximate  subspace  used.  As  further  shown  in  [P5],  (2.4)  is 
sufficient  as  well,  provided  the  remainder  zt  satisfies  “Condition(a)”,  i.e. 

(2.5)  |N|„  <  CO'2 


in  an  appropriate  norm  //.  We  have  derived  several  related  general  theorems  relating 
locking  and  robustness  in  different  norms  in  [P5|.  In  particular,  we  have  shown,  by 
means  of  several  examples,  that  locking  depends  upon  several  factors,  like  the  regularity 
of  the  solution,  the  norm  in  which  the  error  is  calculated,  etc. 

The  general  theory  summarized  above  has  been  used  by  us  to  analyze  several  ex¬ 
amples,  including: 


(1)  In  [P5],  we  have  investigated  the  case  of  heat  flow  through  highly  anisotropic 
materials,  where  t  represents  the  ratio  of  conductivities  in  two  mutually  perpen¬ 
dicular  directions.  “Condition  (a)”  does  not  hold  for  this  problem.  One  of  the 
results  obtained  is  that  locking  cannot  be  avoided  for  the  h  version  unless  the 
mesh  is  properly  aligned,  a  result  similar  to  the  one  for  shells  in  [5].  Also,  it  is 
shown  that  the  p  version  is  robust  for  this  problem. 
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(2)  The  problem  of  locking  in  the  Timoshenko  beam,  which  was  analyzed  in  [1],  [3] 
is  treated  in  [P5]  using  our  theory.  In  particular,  it  is  shown  that  Condition  (o) 
holds,  so  that  (2.4)  is  necessary  and  sufficient.  For  the  h  version,  we  get  locking 
for  any  p  while  for  the  p  version,  no  locking  occurs. 

(3)  In  [P6],  we  derive  various  results  for  Poisson  ratio  locking  (i/  — »  1/2)  when 
straight-sided  elements  in  2-d  are  used.  In  particular,  we  establish  Condition 
(a)  for  this  problem.  The  fact  that  locking  is  avoided  using  the  p  version  or  using 
the  h  version  with  triangular  elements  of  degree  >  4  (obtained  originally  in  [10], 
[7]  respectively)  are  shown  to  follow  immediately  from  our  theory.  Moreover,  it 
is  rigorously  established  that  for  p  <  3,  locking  leads  to  exactly  an  O(h)  loss  of 
convergence  with  triangular  elements.  We  also  establish  various  new  results  for 
two  types  of  rectangular  elements.  Extensions  to  the  3-d  case  are  also  considered. 

As  noted  experimentally  in  [9],  if  the  elements  are  curved,  then  the  effect  of 
locking  is  much  more  severe.  We  have  analyzed  this  case  both  when  the  elements 
can  be  represented  by  polynomial  mappings  and  analytic  mappings.  Our  results 
for  the  p  version  show  that  the  there  is  no  locking  using  polynomial  mappings 
but  the  rate  of  convergence  is  “shifted”,  i.e.  it  is  C(p  —  s)~°  instead  of  Cp~a. 
We  can  explicitly  characterize  s  in  terms  of  the  degree  of  the  mapping.  We  will 
present  these  results  in  a  forthcoming  paper.  A  sample  result  was  presented  in 
[P4],  where  locking  effects  for  the  h  and  p  version  were  compared. 

(4)  We  have  also  considered  the  Reissner-Mindlin  plate  problem.  It  is  known  that 
the  solution  has  a  boundary  layer  for  a  finite  plate  [2],  so  that  Condition  (a) 
may  not  hold  for  H  regular  enough.  We  therefore  consider  the  cases  of  a  plate 
with  periodic  boundary  conditions  and  establish  Condition  (a)  in  this  case.  Once 
again,  we  show  that  the  p  version  is  free  from  locking,  while  the  h  version  shows 
locking  for  any  p  ( 0(h)  loss  of  convergence  for  triangles).  These  and  other  results 
arc  presented  in  [PI 4]. 

(5)  Various  models  for  arches  has  been  discussed  in  [12].  We  have  established  that 
Condition  (a)  holds  in  all  these  cases.  This  enables  us  to  investigate  the  locking  of 
various  standard  finite  element  schemes,  which  will  be  presented  in  a  forthcoming 
paper. 

In  the  above  examples,  we  have  mainly  dealt  with  standard ,  as  opposed  to  mixed 
finite  element  schemes,  since  typically  the  former  arc  the  only  ones  available  in  terms 
of  most  commercial  codes  and  are  consequently  of  great  practical  interest.  Although 
several  kinds  of  locking  problems  have  been  investigated  by  us  above,  there  remain 
important  cases  that  do  not  fit  into  the  general  structure  studied  so  far.  These  will  be 
the  subject  of  future  research. 
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Some  of  the  above  results  on  locking  were  presented  at  an  invited  lecture  at  the 
Numerical  Methods  for  Elliptic  Systems  Workshop,  Helsinki  Institute  of  Technology, 
Espoo,  Finland  (October  25-28,  1989).  Also,  a  minisymposium  entitled  “The  problem 
of  locking  in  finite  element  analysis”  was  organized  by  the  prinicipa!  investigator  as 
part  of  ICIAM  91  (Washington,  D.C.,  July  8-12,  1991),  at  which  various  aspects  of 
locking  were  discussed. 

2.2.  Results  on  the  p  and  h  —  p  versions  of  the  FEM 

The  main  emphasis  of  the  research  done  under  a  previous  grant  (AFOSR  85-0322)  was 
to  establish  various  basic  theoretical  results  for  the  p  and  h  —  p  versions  of  the  FEM. 
These  included  optimal  approximation  results  for  various  types  of  p  and  h-p  subspaces, 
the  numerical  approximation  of  singularities  by  these  methods,  the  treatment  of  in¬ 
homogeneous  Dirichlet  boundary  conditions  by  the  p  version  and  applications  to  the 
boundary  element  method  (BEM)  for  integral  equations.  These  results  were  presented 
at  an  invited  lecture  at  ICOSAIIOM  ’89  at  Villa  Olmo,  Italy  (June  26-29,  1989)  and 
have  been  summarized  in  the  survey  (P2j. 

During  the  current  grant,  some  additional  questions  have  been  investigated.  One  of 
these  is  the  characterization  of  the  L2  error  for  the  p  version,  over  non-convex  domains. 
It  was  shown  in  [11]  that  for  the  h  version,  one  does  not  get  an  O(h)  improvement 
over  the  Hl  error  in  this  case.  This  is  due  to  the  fact  that  the  usual  duality  argument 
breaks  down  when  the  domain  is  non-convex.  In  contrast,  we  have  shown  in  [P8]  that 
the  p  version  gives  the  full  0(p~x)  improvement,  by  using  a  modified  duality  argument. 
Negative  norm  estimates  and  computational  results  (dealing  with  the  sharpness  of  our 
estimates)  arc  also  presented. 

In  [8],  we  investigated  the  stability  and  convergence  properties  of  the  Raviart- 
Thomas  and  Brezzi-Douglas-Marini  spaces  in  the  context  of  the  mixed  p  version  (and 
h-p  version)  FEM.  We  showed  stability  and  optimal  convergence  (up  to  an  c  >  0)  in 
terms  of  p.  In  [P3],  we  have  extended  our  analysis  of  the  Raviart-Thomas  spaces  to  the 
case  of  quasilincar  problems  and  proven  optimal  convergence  once  more  for  the  mixed 
p  version. 

The  p  version  of  the  BEM  for  3-d  problems  has  been  investigated  in  joint  work  by  the 
principal  investigators.  Various  optimal  approximation  results  for  the  different  types 
of  singularities  that  may  occur  in  this  situation  have  been  established.  In  particular, 
we  present  convergence  results  for  vertex,  edge  and  combined  edge- vertex  singularities 
in  the  report  [P9]  (a  more  detailed  paper,  containing  improved  versions  of  these  results 
together  with  their  proofs,  is  in  preparation)': 

Finally,  we  have  been  investigating  the  problem  of  numerical  quadrature  for  the  p 
version.  This  is  one  of  the  last  remaining  basic  problems  for  the  p  version,  and  one  of 
great  practical  importance,  since  for  3-d  codes  (like  STRIPE),  approximately  50%  of 
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the  cost  of  computation  is  from  calculating  the  integrals  in  the  various  matrices.  In  [P7], 
we  have  established  that  the  asymptotic  order  of  convergence  for  the  source  problem 
remains  unchanged  when  commonly  used  types  of  quadrature  are  employed  instead  of 
exact  integration.  We  also  discuss  situations  where  overintegration  may  be  necessary  to 
recover  this  rate  of  convergence,  when  distorted  elements  are  used.  The  paper  contains 
various  computational  results  related  to  the  use  of  mapped  elements.  In  [P 12],  the 
corresponding  results  (including  an  optimal  convergence  rate)  are  established  for  the 
eigenvalue  problem.  In  addition  to  the  results  reported  in  these  references,  various 
additional  results  have  been  obtained  in  collaboration  with  a  Ph.D.  student,  Chang 
Kim,  on  the  effect  of  numerical  integration  on  lower  order  norms. 

To  summarize,  at  this  stage,  most  of  the  basic  theoretical  work  for  the  p  and  h  —  p 
versions  has  been  completed.  This  body  of  work  will  be  reported  in  a  contribution 
being  prepared  (jointly  with  I.  Babuska)  for  the  “Handbook  of  Numerical  Analysis” 
(P.  G.  Ciarlet  and  J.  L.  Lions,  editors).  A  second  survey  [P 13],  which  explains  the  main 
differences  between  the  h  and  the  p  (and  h  —  p)  versions  will  be  submitted  to  SIAM 
Review  shortly. 

2.3._Results  obtained  by  Dr.  C.  Schwab 

The  co-principal  investigator  has  also  completed  research  on  several  other  topics,  por¬ 
tions  of  which  have  been  supported  by  AFOSR  grant  89-0252. 

The.  regularity  of  solutions  to  the  equations  of  3-d  elasticity  and  hydrodynamics 
in  domains  with  edges  and  vertices  is  governed  by  the  singular  exponents  associated 
with  them.  While  the  edge  singularities  are  easily  obtained  from  2-d  theory,  the  vertex 
exponents  are  difficult  to  estimate.  In  [Cl],[C2]  we  proved  sharp  estimates  for  these 
exponents  in  dependence  on  the  geometry  of  the  domain  and  established  a  (nonlinear) 
variational  principle  for  them.  This  implies  in  particular  a  H 2  regularity  result  for 
convex  polyhedra  in  three  dimensions  for  the  Lame-  and  Stokes  equations.  Singular 
exponents  for  rotational  cones  have  been  characterized  by  a  separation  of  variables 
technique,  generalizing  work  by  Ba2ant  and  Keer,  and  we  computed  these  exponents 
numerically  in  [C3] . 

Furthermore,  we  analyzed  and  solved  the  problem  of  numerical  integration  for  the 
singular  surface  integrals  arising  in  boundary  element  methods.  In  [C4],  we  analyzed 
and  classified  the  general  form  of  the  integrands  arising  in  BEM  (our  analysis  covers 
all  operators  that  can  possibly  arise  in  BEM).  Based  on  this,  in  [C5]  we  developed 
numerical  quadrature  schemes  and  proved  error  estimates  in  dependence  on  the  size  of 
the  integration  domain  and  the  number  of  quadrature  points.  In  forthcoming  work  [C7], 
we  will  present  numerical  examples  and  demonstrate  the  use  of  symbolic  manipulation 
in  the  regularization  of  the  hyper  singular,  finite  part  surface  integrals.  Moreover,  we 
also  present  estimates  of  the  quadrature  error  for  the  p-BEM  there. 
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We  also  investgated  in  {C6j  hierarchical  plate  models  for  homogeneous,  isotropic 
plates.  There,  using  the  Fourier  analysis  of  plate  models  developed  in  [6],  we  obtained 
optimal  shear  correction  factors  for  the  Rcissner  Mindlin  and  the  next  higher  model 
(the  (i,l,2)-model)  in  the  hierarchy.  We  also  showed  that  in  all  models  of  higher  order, 
the  introduction  of  a  shear  correction  factor  downgrades  the  asymptotic  accuracy. 

Furthermore,  work  is  in  progress  regarding  a  coercive  formulation  for  exterior  three 
dimensional  solid  -  fluid  interaction  problems  that  can  then  be  combined  with  dimen¬ 
sional  reduction  techniques  to  obtain  well  posed  boundary  value  problems  for  shells 
immersed  in  a  fluid. 
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